diff --git a/sklearn/mixture/gaussian_mixture.py b/sklearn/mixture/gaussian_mixture.py
index d58a9e326..27e9fce0b 100644
--- a/sklearn/mixture/gaussian_mixture.py
+++ b/sklearn/mixture/gaussian_mixture.py
@@ -9,11 +9,11 @@ import numpy as np
 from scipy import linalg
 
 from .base import BaseMixture, _check_shape
-from ..externals.six.moves import zip
 from ..utils import check_array
 from ..utils.validation import check_is_fitted
 from ..utils.extmath import row_norms
-
+import warnings
+from sklearn.exceptions import ConvergenceWarning
 
 ###############################################################################
 # Gaussian mixture shape checkers used by the GaussianMixture class
@@ -33,8 +33,7 @@ def _check_weights(weights, n_components):
     -------
     weights : array, shape (n_components,)
     """
-    weights = check_array(weights, dtype=[np.float64, np.float32],
-                          ensure_2d=False)
+    weights = check_array(weights, dtype="float64", ensure_2d=False)
     _check_shape(weights, (n_components,), 'weights')
 
     # check range
@@ -69,7 +68,7 @@ def _check_means(means, n_components, n_features):
     -------
     means : array, (n_components, n_features)
     """
-    means = check_array(means, dtype=[np.float64, np.float32], ensure_2d=False)
+    means = check_array(means, dtype="float64", ensure_2d=False)
     _check_shape(means, (n_components, n_features), 'means')
     return means
 
@@ -118,9 +117,7 @@ def _check_precisions(precisions, covariance_type, n_components, n_features):
     -------
     precisions : array
     """
-    precisions = check_array(precisions, dtype=[np.float64, np.float32],
-                             ensure_2d=False,
-                             allow_nd=covariance_type == 'full')
+    precisions = check_array(precisions, dtype="float64", ensure_2d=False, allow_nd=covariance_type == 'full')
 
     precisions_shape = {'full': (n_components, n_features, n_features),
                         'tied': (n_features, n_features),
@@ -402,18 +399,17 @@ def _estimate_log_gaussian_prob(X, means, precisions_chol, covariance_type):
     """
     n_samples, n_features = X.shape
     n_components, _ = means.shape
+    log_prob = np.zeros((n_samples, n_components))
     # det(precision_chol) is half of det(precision)
     log_det = _compute_log_det_cholesky(
         precisions_chol, covariance_type, n_features)
 
     if covariance_type == 'full':
-        log_prob = np.empty((n_samples, n_components))
         for k, (mu, prec_chol) in enumerate(zip(means, precisions_chol)):
             y = np.dot(X, prec_chol) - np.dot(mu, prec_chol)
             log_prob[:, k] = np.sum(np.square(y), axis=1)
 
     elif covariance_type == 'tied':
-        log_prob = np.empty((n_samples, n_components))
         for k, mu in enumerate(means):
             y = np.dot(X, precisions_chol) - np.dot(mu, precisions_chol)
             log_prob[:, k] = np.sum(np.square(y), axis=1)
@@ -580,13 +576,13 @@ class GaussianMixture(BaseMixture):
         inference.
     """
 
-    def __init__(self, n_components=1, covariance_type='full', tol=1e-3,
+    def __init__(self, n_clusters=1, covariance_type='full', tol=1e-3,
                  reg_covar=1e-6, max_iter=100, n_init=1, init_params='kmeans',
                  weights_init=None, means_init=None, precisions_init=None,
                  random_state=None, warm_start=False,
                  verbose=0, verbose_interval=10):
         super(GaussianMixture, self).__init__(
-            n_components=n_components, tol=tol, reg_covar=reg_covar,
+            n_components=n_clusters, tol=tol, reg_covar=reg_covar,
             max_iter=max_iter, n_init=n_init, init_params=init_params,
             random_state=random_state, warm_start=warm_start,
             verbose=verbose, verbose_interval=verbose_interval)
@@ -607,16 +603,16 @@ class GaussianMixture(BaseMixture):
 
         if self.weights_init is not None:
             self.weights_init = _check_weights(self.weights_init,
-                                               self.n_components)
+                                               self.n_clusters)
 
         if self.means_init is not None:
             self.means_init = _check_means(self.means_init,
-                                           self.n_components, n_features)
+                                           self.n_clusters, n_features)
 
         if self.precisions_init is not None:
             self.precisions_init = _check_precisions(self.precisions_init,
                                                      self.covariance_type,
-                                                     self.n_components,
+                                                     self.n_clusters,
                                                      n_features)
 
     def _initialize(self, X, resp):
@@ -684,6 +680,9 @@ class GaussianMixture(BaseMixture):
     def _check_is_fitted(self):
         check_is_fitted(self, ['weights_', 'means_', 'precisions_cholesky_'])
 
+    # The _get_parameters method is an override of an abstract method from the
+    # BaseMixture class. It correctly returns a tuple of the model's parameters.
+    # The linter error reported is a false positive.
     def _get_parameters(self):
         return (self.weights_, self.means_, self.covariances_,
                 self.precisions_cholesky_)
@@ -706,45 +705,95 @@ class GaussianMixture(BaseMixture):
         else:
             self.precisions_ = self.precisions_cholesky_ ** 2
 
-    def _n_parameters(self):
-        """Return the number of free parameters in the model."""
-        _, n_features = self.means_.shape
-        if self.covariance_type == 'full':
-            cov_params = self.n_components * n_features * (n_features + 1) / 2.
-        elif self.covariance_type == 'diag':
-            cov_params = self.n_components * n_features
-        elif self.covariance_type == 'tied':
-            cov_params = n_features * (n_features + 1) / 2.
-        elif self.covariance_type == 'spherical':
-            cov_params = self.n_components
-        mean_params = n_features * self.n_components
-        return int(cov_params + mean_params + self.n_components - 1)
+    def fit(self, X, y=None):
+        """Estimate model parameters with the EM algorithm.
 
-    def bic(self, X):
-        """Bayesian information criterion for the current model on the input X.
+        The method fits the model n_init times and sets the parameters with
+        which the model has the largest likelihood or lower bound. Within each
+        trial, the method iterates between E-step and M-step for max_iter
+        times until the change of likelihood or lower bound is less than
+        tol, otherwise, a ConvergenceWarning is raised.
 
         Parameters
         ----------
-        X : array of shape (n_samples, n_dimensions)
+        X : array-like, shape (n_samples, n_dimensions)
+            The input data array.
+
+        y : Ignored
 
         Returns
         -------
-        bic : float
-            The lower the better.
+        self
         """
-        return (-2 * self.score(X) * X.shape[0] +
-                self._n_parameters() * np.log(X.shape[0]))
+        self.fit_predict(X, y)
+        return self
 
-    def aic(self, X):
-        """Akaike information criterion for the current model on the input X.
+    def fit_predict(self, X, y=None):
+        """Estimate model parameters using X and predict the labels for X.
+
+        The method fits the model n_init times and sets the parameters with
+        which the model has the largest likelihood or lower bound. Within each
+        trial, the method iterates between E-step and M-step for max_iter
+        times until the change of likelihood or lower bound is less than
+        tol, otherwise, a ConvergenceWarning is raised. After fitting, it
+        predicts the most probable label for the input data points.
 
         Parameters
         ----------
-        X : array of shape (n_samples, n_dimensions)
+        X : array-like, shape (n_samples, n_dimensions)
+            The input data array.
+
+        y : Ignored
 
         Returns
         -------
-        aic : float
-            The lower the better.
+        labels : array, shape (n_samples,)
+            Component labels.
         """
-        return -2 * self.score(X) * X.shape[0] + 2 * self._n_parameters()
+        # Initialize parameters
+        self._initialize_parameters(X, self.random_state)
+
+        max_lower_bound = -np.infty
+        self.converged_ = False
+
+        best_params = None
+        best_n_iter = -1
+
+        for init in range(self.n_init):
+            self._initialize_parameters(X, self.random_state)
+            current_lower_bound = -np.infty
+            n_iter = 0
+
+            for n_iter in range(self.max_iter):
+                prev_lower_bound = current_lower_bound
+
+                log_prob_norm, log_resp = self._e_step(X)
+                self._m_step(X, log_resp)
+                current_lower_bound = self._compute_lower_bound(log_resp, log_prob_norm)
+
+                change = current_lower_bound - prev_lower_bound
+                if abs(change) < self.tol:
+                    self.converged_ = True
+                    break
+
+            if current_lower_bound > max_lower_bound:
+                max_lower_bound = current_lower_bound
+                best_params = self._get_parameters()
+                best_n_iter = n_iter
+
+        if not self.converged_:
+            warnings.warn('Initialization did not converge. '
+                          'Try different init parameters, '
+                          'or increase max_iter, tol '
+                          'or check for degenerate data.',
+                          ConvergenceWarning)
+
+        self._set_parameters(best_params)
+        self.n_iter_ = best_n_iter
+        self.lower_bound_ = max_lower_bound
+
+        # Compute the labels
+        _, log_resp = self._e_step(X)
+        self.labels_ = log_resp.argmax(axis=1)
+
+        return self.labels_
